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Abstract: We give a preliminary report on one of the tests we have performed of 
a fuU non-axisymmetric general relativistic code. The test considered here concerns 
the numerical evolution of vacuum non-axisymmetric gravitational waves and their 
comparison at low amplitudes with theoretical waveforms obtained from linearised 
theory. 
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I. Introduction 

We are currently carrying out a number of detailed tests of a full non-axisymmetric 
general relativistic code which solves the coupled Einstein and hydrodynamic equations 
using the '3-M' ADM formalism in the radial gauge with mixed slicing [1]. The code 
has been written in particular to obtain directly the gravitational wave emission from 
non-axisymmetric rotating gravitational collapse. The test we report here concerns the 
numerical evolution of vacuum non-axisymmetric gravitational waves and their comparison 
at low amplitudes with theoretical waveforms obtained using linearised theory. This report 
should be considered preliminary, since the code is still under active development. 

We consider the propagation of non-axisymmetric I = 2 and / = 4, and different 
m, linearised gravitational waves. We are able to use essentially analytic predictions for 
these waves with which to compare to the numerical solutions. To obtain these, we take 
an even mode linearised axisymmetric solutions for a given 1 = 1^, and rotate this by 7r/2 



E-mail address: richard@fiz.uni-lj.si and norver@starkl.freeserve.co.uk 



1 



so that the symmetry axis now hes in the equator. This generates a non-axisymmetric 
solution which is a hnear combination of all I up to lyj and all even m. (This solution also 
satisfies the symmetry and regularity conditions we have assumed). Under this rotation, 
both even and odd waves are generated. The solution is determined by an arbitrary wave 
amplitude F{t—r) or F{t-\-r) for the retarded or advanced solutions respectively. Although 
each solution is by itself singular at the origin, the combination of (advanced — retarded) 
solution is non-singular everywhere. We choose a particular function F; construct initial 
data at t = 0; and evolve this initial data and observe the outgoing waveform leaving the 
outer grid radius. We then compare this with the analytic solution. 



II. Radial gauge and mixed slicing 

We will not list here the full non-axisymmetric '3+1' equations for our gauge choice. 
These can be found in ref. [2] and will be more fully discussed in a forthcoming paper [1] . 
We give here only the briefest details in order to define the variables used in our wave 
solution. 

We use spherical spatial coordinates (r,^,^). The spatial metric hij = 1,2,3) 
in the radial gauge [1-4] is determined from the conditions that the off-diagonal (r,6') (r,(/») 
components of the metric be zero, and that the angular spatial determinant be fixed to its 
Euclidean value. Thus we have: 

hre = 0, (la) 
hr4> = 0, (lb) 

{heeh^<l> - h^^) = r'^ sin^ 0. (Ic) 

The 3-dimensional metric is then determined by the three functions A, rj, and ^, each 
functions of (t, r, 9, 0), and the full 4-dimensional metric has the form: 

ds^ = -{N^ - N'Ni)dt^ - 2Nidx'dt + A^dr^ + r^B'^dO^ + r^B^{smed(l) + ^dOf, (2) 

with = 1 + rj, and where rj and ^ are respectively the even and odd dynamical degrees 
of freedom in the sense that they tend (for this linearised test, and within a sign) at large 
radii to the even and odd transverse traceless amplitudes h+,hx- N is the lapse; N'^ the 
shift vector raised and lowered using hij. 

The components of the extrinsic curvature we use are based on the orthonormal 
components = 1,2,3) which are the components of Kij projected onto the 

orthonormal triad of basis vectors: 

el) = [A-\0,0]; e\^^ = [0, B/r, -^B/{rsm9)]; ejg) = [0, 0, l/(5rsin^)]. (3) 
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The computed geometrical variables are (we do not discuss here the accompanying 
hydrodynamic variables) : 

Spatial metric : A; rj; ^; (4a) 
Extrinsic curvature: i^i = = 5i^(i)(2)/ sin^; 

= 5K(i)(3)/sin^; K+ = ^B^K^^^^^^ - K^2)(2)); = B^K^2)(3); (46) 

Lapse and shift : N; = 7V7r; G = smO; iV^. (4c) 

Following ref. [4] we use a mixed foliation condition consisting of a linear combination 
(depending on radius) of polar and maximal slicing: we set 

tTace{Ki^) = (1 - C(r))K/, (5) 

where < C(r) < 1 is a chosen smooth function of radius such that C(0) = 1 and 
C(r > ro) = (corresponding to maximal slicing at the origin and polar slicing for radii 
larger than a chosen transition radius ro). The actual form of C(r) used is: 

C(r) = [1 - {r/rofr, (6) 

with ro = ro{t), n chosen values. 

Regularity at the origin is determined by requiring the components in the 'cartesian' 
{x,y,z) system of coordinates be expandable in non-negative integer powers of {x,y,z) 
(consistent with the symmetry assumed) [4]. {x,y,z are quasi-cartesian coordinates related 
to r,^,0 in the usual manner). This gives, in the linearised case (with no mean origin 
azimuthal shift rotation), the following values for the origin shift and extrinsic curvature: 

/3'-(r = 0) = {/c+(3cos2 6'-l) + (/c_cos20-/cxsin20)sin2 6'}, (7a) 

G{r = 0) = {-3A;+ - sin 20 + A;_ cos 20} cos 9, (76) 

7V<^(r = 0) - {-k- sin 20 - fcx cos 20}, (7c) 
3 1 

K+{r = 0) = sin^ 6* + -(1 + cos^ 0)[-k^ sin 20 + /c_ cos 20], (8a) 

(r = 0) = cos d[k- sin 20 + fcx cos 20] , (86) 

Ki{r = 0) = /c+(l - 3cos^0) + sin2^[-/c_ cos20 + fcx sin 20], (8c) 

K2{r = 0) = cos6'[3/c+ - /c_ cos 20 + /cx sin 20], (8d) 

iiTg (r = 0) = [k- sin 20 + /cx cos 20] . (8e) 
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where k+{t), k-{t), kx{t) are time dependent coefficients. 
Also we have near the origin: 

V = 0{r% (9) 
N = N(r = 0) + O{r^). (10) 

III. Axisymmetric linearised gravitational waves 

The axisymmetric solution we rotate is based on one found by Bardeen [5,6] (see 
also ref.[7] where axisymmetric comparison tests are also shown). Bardeen's solution cor- 
responds to a general I even mode Teukolsky wave [8] transformed from the transverse- 
traceless gauge to the radial gauge with maximal slicing. For completeness, we have 
generalised Bardeen's solution (which is for maximal slicing) to the case of mixed slicing 
by obtaining the linearised axisymmetric equations in the radial gauge with mixed slicing. 
As we explain below however, the original Bardeen solution is actually sufficient, even for 
mixed slicing, for the comparison with the computed outgoing waves. 

Following Bardeen [5,6], we write the even linearised solution for each I >2 in the 

form: 



(A - 1) = aiPio, 


(11a) 




(116) 




(11c) 


Ki = -diPio, 


(lid) 


K2 = —-^biPio,a; , 


(lie) 


i^3 = 0. 


(11/) 


K+ = -^q(1 - x'^)Pio,xx , 


(11^) 


=0, 


(ll/i) 




(llz) 


N'l' = 0, 


(Hi) 


= M + ^^Kl - C(r))]Pzo, 


(Ilk) 


N = niPio, 


(110 
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where • = (^)r- In this form, the evolution equation for rj [1,2] is automatically satisfied. 
The equation for P'^ [1,2] is also satisfied. 

The momentum constraints for Ki and K2 [1,2] require: 

bl = -I(7^^"'[^'«^]- -[^(1 - C(^))M,r }, (12) 

= (/_i|(/ + 2) ^^"'^^'^'J'^ ^""^ + " (13) 
The shift equation for G [1,2] provides an ordinary differential equation for ff. 

2r/i„ +/(/ + l)fi = -rai,r -6ai + 2[r(l - C(r))az]„ + 1)(1 - C(r))ai. (14) 

The origin boundary conditions (7), (8) (in the axisymmetric limit) and (9), give for the 
origin boundary condition for ff. 

;,(. = 0) = {--(^ = °)' (15) 

The evolution equation for A [1,2] gives: 

1 

ai = ai + fi + rfi,r -^li'^ - C{r))rai],r , (16) 
or equivalently using (14) for ff. 

ai = -2ai - ^rafr -^{1 - 1)(/ + 2)fi + ^[(1 - C{r))rai],r + 1)(1 - C{r))ai. (17) 

The evolution equation for Ki [1,2] gives (using (16), (17)): 

r'^di = r'^Ufrr +r'^ai,rr +Qrai,r -{l + 3){l- 2)ai - [r^(l - C{r))ai],rr , (18) 
while the lapse equation [1,2], using (18), gives: 

r\i,rr +2rni,r -l{l + l)ni = (1 - C{r))r'^au (19) 
with origin boundary conditions (10): 

n;(r = 0) = 0; n^,^ (r = 0) = 0. (20) 

The remaining equations: the Hamiltonian constraint [1,2], and the K2 evolution equation 
[1,2], are then also satisfied by these solutions. 
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For the case of maximal slicing (C(r) = 1) we can explicitly solve (18), (19) for ai, 
ni and recover Bardeen's solution. In this case, the solution for ni is simply: 



ni = 0. (21) 

The solution for ai for the case of a retarded solution (where t occurs only in the form of 
retarded time (t — r)) is found using the form: 



= Y.CnT-^^^'^Fl'-^\t-T), (22) 



where indicates the fc'th derivative with respect to the arbitrary generating function 
Fi for the Tth mode. Substitution into (18) gives for the coefficients c„: 



_ (Z + n)(/-n+l) 



leading to the retarded solution: 



I 



= E 7,^l^.r-^''^'^Fl'-\t - r). (24) 
^2«n!(Z-n)! i \ ) ^ > 

n=0 ^ 



The advanced solution is: 

I 

2'^n\{l-n)\ 



n=0 



while the non-singular (advanced — retarded) solution is: 



n=0 



Substituting the non-singular solution (26) for ai into (14) (with C(r) = 1) we find 

for //: 

-r[{-l)''F[^-''+^\t+r)+F[^~''+^\t-r)] + {n-3)[{-l)''Fl^^^^^ (27) 

The integral can be evaluated either directly, or by reducing it using integration by parts 
to recurrence relations involving the Fl^\t ± r). 
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The outgoing waveform at asymptotically large radii of the transverse traceless 
amplitudes is: 

''+ = -(I3TwTT)(i + 2)'-"'^''"'('-'-«i--')^'°- P8) 

/ix = 0. (29) 



IV. Generation of non-axisymmetric linearised gravitational waves 

We now rotate this solution, with the symmetry axis of the wave rotated into 
the equator, to generate a non-axisymmetric solution. The relation between {0,4)) to the 
rotated {0',(f)') coordinates aligned with the symmetry axis of the wave (which when rotated 
lies in the equator at an angle a to the x axis of the {9,(f)) system, corresponding to x' = —z; 
y' = y cos a — X sin a; z' = x cos a + y sin a; r' = r) is: 

sin 9' cos 4>' — — cos 9, (30a) 

sin 9' sin </>' = sin 9 sin(</) — a), (306) 
cos6'' = sin6'cos(0 — a). (30c) 
Under this rotation, to linear order, the geometrical variables become: 

r]{t,r,9,(t>)= I' [cos 2(0 - a) - - sin^ 6>(1 + cos 2(0 - a))], (31a) 
sin 9' 2 

e(t, r, 9, 0) = "^'^^^l^p cos 9 sin 2(0 - «), (316) 
sm 9' 

G{t, r, 9, 0) = cos ^ J(l + cos 2(0 - a)), (31c) 

COS u 

N'^{t, r, 9, 0) = ^^M^ 1 sin 2(0 - a), (31d) 

COS (/ ^ 

(i^{t,r,9,<\>) = r{t.r,9'), (31e) 

A(t,r,^,0) = A'(t,r,^'), (31/) 

iV(t,r,^,0) = iV'(t,r,^'), (31^) 
K' (t r 9') 1 

K+{t,r,9,(f)) = +^ ; ' ^ -[-sin^6> + (l + cos^6>)cos2(0-a)], (31/i) 
sin ^' 2 

K' it r 9') 

{t, r, 9, 0) = + 3 ' ^ cos ^ sin 2(0 - a), (31i) 
sin 9' 
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K2{t, r, e, 4>) = „7 ^'^ I COS ^(1 + cos 2(0 - a)), {31k) 

COS u Zi 

m, r, 9, cf>) = ^f^jf^ \ sin 2(0 - a). (31/) 

Applying these rotated values to the non-singular solution (26), (27) and expanding 
around the origin, we find the following values for the origin coefficients k+{t), k-{t), kx {t) 
in (7), (8) when / = 2: 

k+ = -^Fitl{t)\r=o, (32a) 



A;- = ^FS(t)|r=oCos2a, (326) 
fcx = -^F^^l{t)\r=osm2a. (32c) 



For I > 2 these coefficients are zero. 

The I = 2 solution after rotation corresponds to the following linear combination of 

waves: 

r7 = ^2(r,t)[-^P22+(3- ^P22)cos2(0-a)], (33a) 

C = ?72(r,t)3Piosin2(0-a), (336) 
The I = 4 rotated solution corresponds to the linear combination: 

3 15 5 1 5 1 

^ = ^4(r,t){-P42 + (-Y + 2^'22+2A2)cos2(0-a) + (-P22 + gP42)cos4(0-a)}, (34a) 

15 7 7 

e = 774 (r, t){{-—Pio + ^P32) sin 2(0 - a) + -P32 sin 4(0 - a)}, (346) 

where in the above: 

^Kr,*) = [Q + ^|^]- (35) 



Expanding Fi near the origin, we find the leading order non-zero terms of "fji for 
/ = 2,4 to be: 

m = -j^r'F[ll{t)\r=o (r^O), (36) 

^4 = -^r^FS(t)|.=o (r-0). (37) 

This leading order behaviour is consistent with our assumption of cartesian integer ex- 
pandability near the origin. 
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The outgoing form for fji{r,t) at asymptotically large radii is: 

Substitution of (38) into (33), (34) gives the outgoing waveform at asymptotically large 
radii. The transverse traceless amplitudes (to linear order) are then simply: 

h+ = -ri, (39) 

/ix = e (40) 



V. The numerically evolved solution 

For the / = 2 wave test, we choose the wave amplitude function Fi to be: 

i^(«=2)(y) = ayexp(-/). (41) 

For / = 4 we choose: 

F(z=4)(y) =ay3exp(-y2), (42) 

where a is the (small) amplitude of the wave. The extra factor of in (42) is to simulate 
more complex radial waveform structure. 

The moments F^'^'^ in the solutions (26), (27) for our choices (39), (40) for Fi can be 
automatically generated essentially analytically by repeated application of the recurrence 
relation: 

yPexpi-y') ^ [pyP-' - 2yP+') exp{-y^), (43) 

for the mapping on differentiation of the coefficients multiplying the powers of y. Hence 
if required, the whole linearised solution (for maximal slicing) can be obtained essentially 
'analytically' throughout the full spacetime for comparison with numerical code solutions. 

An important point to note is that although our 'analytic' solution is for maximal 
slicing, we can still use it for comparison with the outgoing solution at large radii obtained 
numerically with mixed slicing. This is because firstly, the choice (41), (42) corresponds 
at t = to time symmetric initial data with {Ki,K2,K2,K^,Kx) all zero. Hence initially 
the mixed and maximal slicing are equivalent. Secondly during the evolution at large r 
where we monitor the outgoing wave, h+,hx are gauge invariant while the coordinates 
times agree (differing only by terms ~ a^/r). 

h+ and hx can be spectrally represented (to linear order) in the following manner 

[1]: 

h+ = (-E'axis cos(20) + i:)axissin(2(;i;))) 
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- ^ {ai + aio cos(20) + bio sin(20))P/2 

1=2,4,6... 

- (a/mCos((m + 2)0) + 6/^sin((m + 2)0))P/^, (44a) 

i=2,4,6,... 
m=2,4,...,l 

hy, = (Daxis C0S(2^) +£;axisSin(2(/)))Pi 
1=3,5,7... 

+ Y WmCOs((m + 2)0) + 6;^sin((m + 2)0))Pi^, (446) 

1=3,5,7... 
m=2,4,..,l-l 

where the spectral amplitudes -Daxis^-^'axis! O'l] O'loMo'i o-lmMm ctnd their primed counter- 
parts, arc functions of radius and time. These expansions are used to spectrally resolve 
the angular dependence of the numerically obtained outgoing waves on the angular grid 
points at each radius. 

For the numerical evolution, we have used an outer grid radius of 10 units; a grid 
of (100x11x16) (r, 0, 0) grid points; and an amplitude a = 10~^ for / = 2 and a = 10~^ for 
1 = 4. vq is 1.4534 units, and n = 2.5 (see eq.(6)). We have taken a = 0. (Details of the 
numerical code will be given in a separate paper [9]). The expected non-zero asymptotic 
amplitudes (for these values) are, for I = 2: 

[£^axis,a2,a2o] = ^ [6,-1,-1] r-iF((£2)(^-0 (r ^ oo), (45a) 

where: 

4=2)^^) ^ aexp(-y2)(60y - 80y^ + 16^^), (456) 

while for i = 4: 

[£;axis, ^4, ^20, ^40, ^22, ^42, 630, 633] = [-60, 3, 20, 4, 20, 1, 14, 7] 

r-^F^il^^{t-r) (r^oo) (46a) 

with: 

^(£4)(^) = aexp(-|/2)(2520|/ - 7560y^ + 5040y^ - 1056y^ + 64y^). (466) 

Figures 1-11 show a comparison of the numerically computed non-zero spectral 
outgoing amplitudes versus the theoretically expected asymptotic waveform amplitudes 
(45), (46). Figures 1-3 show results for the / = 2 test; figures 4-11 results for the / = 4 test. 
The numerical waveforms are measured at the outer grid radius of 10 units. Reasonable 
agreement is obtained, although the later part of the numerical waveforms tend to damp 
too slowly and eventually show irregularity. The present preliminary setup of our code is 
unable to run much further than shown because of stability difficulties. Further work is 
currently underway to overcome these problems. 
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The tests presented here are derived from rotated axisymmetry solutions, and are 
not the most general non-axisymmetric test possible. More generally, radial gauge coordi- 
nates satisfying the origin regularity conditions will not have all metric and shift variables 
which explicitly fall-off asymptotically at large radii [1] . These more general types of tests 
will be dealt with in a separate paper [9] . 
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Figure captions: 



Fig.l: I = 2 wave test; -Eaxis amplitude (see eqs.(44),(45)) vs. retarded time (t — r); numer- 
ical and theoretical curves. Comparison of the measured outgoing amplitude (full curve) 
from the full numerical code with the theoretical linearised prediction (broken curve). The 
amplitude is normalised by the small initial amplitude a. Show is the variation of the 
amplitude vs. retarded time. The numerical solution is obtained at the outer grid radius 
of 10 units. The theoretical curve is that for asymptotically large radii. 

Fig. 2: / = 2 wave test; a2 amplitude (see eqs.(44),(45)) vs. retarded time (t — r); numerical 
and theoretical curves, (see also Fig.l caption). 

Fig. 3: / = 2 wave test; 020 amplitude (see eqs.(44),(45)) vs. retarded time (t—r); numerical 
and theoretical curves, (see also Fig.l caption). 

Fig. 4: / = 4 wave test; i?axis amplitude (see eqs.(44),(46)) vs. retarded time {t — r); 
numerical and theoretical curves, (see also Fig.l caption). 

Fig. 5: I = 4 wave test; 04 amplitude (see eqs.(44),(46)) vs. retarded time (t — r); numerical 
and theoretical curves, (see also Fig.l caption). 

Fig. 6: / = 4 wave test; 020 amplitude (see eqs.(44),(46)) vs. retarded time (t—r); numerical 
and theoretical curves, (see also Fig.l caption). 

Fig.7: / = 4 wave test; a4o amplitude (see eqs.(44),(46)) vs. retarded time (t—r); numerical 
and theoretical curves, (see also Fig.l caption). 

Fig. 8: / = 4 wave test; 022 amplitude (see eqs.(44),(46)) vs. retarded time (t—r); numerical 
and theoretical curves, (see also Fig.l caption). 

Fig. 9: / = 4 wave test; 042 amplitude (see eqs.(44),(46)) vs. retarded time (t—r); numerical 
and theoretical curves, (see also Fig.l caption). 

Fig. 10: / = 4 wave test; 639 amplitude (see eqs.(44),(46)) vs. retarded time {t — r); 
numerical and theoretical curves, (see also Fig.l caption). 

Fig.ll: I = 4 wave test; 632 amplitude (see eqs.(44),(46)) vs. retarded time {t — r); 
numerical and theoretical curves, (see also Fig.l caption). 
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